/**
  SUNNY project, Anyang Normal University, IMP-CAS
  \class TACGCoe
  \brief calculate Clebsch-Gordan coefficient <j1m1j2m2|jm>
  \author SUN Yazhou, asia.rabbit@163.com
  \since 2021/02/28
  \date Last modified: 2021/02/28 by SUN Yazhou
  \copyright 2020-2021 SUN Yazhou
  \copyright CNOK project, Anyang Normal University, IMP-CAS
*/

#include <cmath>
#include <algorithm>
#include "config.h"
#ifdef CATCH_CNOK
#include "catch.h"
#endif
#include "TACGCoe.h"
#include "TAException.h"
#include "TAMath.h"

static auto fn = TAMath::factln; // fn(n) = ln(n!)
static auto mi = TAMath::minus; // mi(n) = (-)^n
using std::min;
using std::max;

// 21 = 0.5, m2 = 0.5
double TACGCoe::CGHalfP(double j1, double m1, double j, double m){
  if(j == j1+0.5) return sqrt((j1+m+0.5)/(2.*j1+1.));
  if(j == j1-0.5) return -sqrt((j1-m+0.5)/(2.*j1+1.));
  // TAException::Error("TACGCoe", "CGHalfP: illegal j=%d.", j);
  return 0.;
} // end member function CGHalfP

// j2 = 0.5, m2 = -0.5
double TACGCoe::CGHalfM(double j1, double m1, double j, double m){
  if(j == j1+0.5) return sqrt((j1-m+0.5)/(2.*j1+1.));
  if(j == j1-0.5) return sqrt((j1+m+0.5)/(2.*j1+1.));
  // TAException::Error("TACGCoe", "CGHalfM: illegal j=%d.", j);
  return 0.;
} // end member function CGHalfM

// j2 = 1, m2 = 1
double TACGCoe::CG1P(double j1, double m1, double j, double m){
  if(j == j1+1) return sqrt((j1+m)*(j1+m+1.)/((2.*j1+1.)*(2.*j1+2.)));
  if(j == j1) return -sqrt((j1+m)*(j1-m+1.)/(2.*j1*(j1+1.)));
  if(j == j1-1) return sqrt((j1-m)*(j1-m+1.)/(2.*j1*(2.*j1+1.)));
  // TAException::Error("TACGCoe", "CG1P: illegal j=%d.", j);
  return 0.;
} // end member function CG1P

// j2 = 1, m2 = 0
double TACGCoe::CG1Z(double j1, double m1, double j, double m){
  if(j == j1+1) return sqrt((j1-m+1.)*(j1+m+1.)/((2.*j1+1.)*(j1+1.)));
  if(j == j1) return m/sqrt(j1*(j1+1.));
  if(j == j1-1) return sqrt((j1-m)*(j1+m)/(j1*(2.*j1+1.)));
  // TAException::Error("TACGCoe", "CG1Z: illegal j=%d.", j);
  return 0.;
} // end member function CG1Z

// j2 = 1, m2 = -1
double TACGCoe::CG1M(double j1, double m1, double j, double m){
  if(j == j1+1) return sqrt((j1-m)*(j1-m+1.)/((2.*j1+1.)*(2.*j1+2.)));
  if(j == j1) return -sqrt((j1-m)*(j1+m+1.)/(2.*j1*(j1+1.)));
  if(j == j1-1) return sqrt((j1+m)*(j1+m+1.)/(2.*j1*(2.*j1+1.)));
  // TAException::Error("TACGCoe", "CG1M: illegal j=%d.", j);
  return 0.;
} // end member function CG1M


inline bool valid(double j, double m){
  if(j < 0. || fabs(m) > j) return false;
  const double d = j - floor(j), dm = m - floor(m);
  if((d != 0. && d != 0.5) || d != dm) return false;
  return true;
}
double TACGCoe::CG(double j1, double m1, double j2, double m2, double j, double m){
  if(!valid(j1,m1) || !valid(j2,m2) || !valid(j,m))
    TAException::Error("TACGCoe", "CG: Invalid input.");
  if(j > j1+j2 || j < fabs(j1-j2) || m != m1+m2) return 0.;

  // accommodate for j1 or j2 = 1 or 0.5.
  if(0.5 == j2) return m2 > 0 ? CGHalfP(j1,m1,j,m) : CGHalfM(j1,m1,j,m);
  if(0.5 == j1) return mi(j1+j2-j)*(m1 > 0 ? CGHalfP(j2,m2,j,m) : CGHalfM(j2,m2,j,m));
  if(1. == j2) return m2>0?CG1P(j1,m1,j,m):(m2?CG1M(j1,m1,j,m):CG1Z(j1,m1,j,m));
  if(1. == j1) return mi(j1+j2-j)*(m1>0?CG1P(j2,m2,j,m):(m1?CG1M(j2,m2,j,m):CG1Z(j2,m2,j,m)));

  // Racah's general formula for CG //
  const double c = fn(j1+j2-j)+fn(j2+j-j1)+fn(j+j1-j2)+
    fn(j1+m1)+fn(j1-m1)+fn(j2+m2)+fn(j2-m2)+fn(j+m)+fn(j-m)-fn(j1+j2+j+1.);
  const double vmin = max({0.,j1+m2-j,j2-m1-j}), vmax = min({j1+j2-j,j1-m1,j2+m2});
  double s = 0.;
  for(double v = vmin; v <= vmax; v++)
    s += mi(v)*exp(-(fn(v)+fn(j1+j2-j-v)+fn(j1-m1-v)+fn(j2+m2-v)+fn(j-j1-m2+v)+fn(j-j2+m1+v)));
  return sqrt(2.*j+1.) * exp(c/2.) * s;
} // end member function CG

double TACGCoe::ThreeJSymbol(double j1, double m1, double j2, double m2, double j, double m){
  return mi(j1-j2-m) * CG(j1,m1,j2,m2,j,-m) / sqrt(2.*j+1.);
} // end member function ThreeJSymbol

#ifdef CATCH_CNOK
TEST_CASE("Clebsch-Gordan Coefficient", "[cg]"){
  CHECK(TACGCoe::CG(5,-2, 7,-5, 10,-7) == Approx(-0.1003426623405196).epsilon(1e-10));
  CHECK(TACGCoe::CG(4,4, 12,3, 10,7) == Approx(0.07541506695341409).epsilon(1e-10));
  CHECK(TACGCoe::CG(0.5,0.5, 12,3, 12.5,3.5) == Approx(0.8).epsilon(1e-10));
  CHECK(TACGCoe::CG(18,13, 37,15, 50,28) == Approx(-0.2799711347740473).epsilon(1e-10));
  CHECK(double(TACGCoe::ThreeJSymbol(3,1, 4,2, 5,-3)-TACGCoe::ThreeJSymbol(4,2, 3,1, 5,-3)+1.)
    == Approx(1.).epsilon(1.e-10));
  CHECK(TACGCoe::CG(0.5,-0.5,0.5,-0.5,1.,-1.) == 1.);
  CHECK(TACGCoe::CG(0.5,0.5,0.5,0.5,1.,1.) == 1.);
  CHECK(TACGCoe::CG(0.5,0.5,0.5,-0.5,1.,0.) == TAMath::Sqrt2()/2.);
  CHECK(TACGCoe::CG(0.5,-0.5,0.5,0.5,1.,0.) == TAMath::Sqrt2()/2.);
  CHECK(TACGCoe::CG(0.5,0.5,0.5,-0.5,0.,0.) == TAMath::Sqrt2()/2.);
  CHECK(TACGCoe::CG(0.5,-0.5,0.5,0.5,0.,0.) == -TAMath::Sqrt2()/2.);

  CHECK(TACGCoe::CG(0.5,0.5,1.,0.,1.5,0.5) == Approx(sqrt(6.)/3.).epsilon(1e-10));
  CHECK(TACGCoe::CG(0.5,0.5,1.,1.,1.5,1.5) == Approx(1.).epsilon(1e-10));
  CHECK(TACGCoe::CG(3.,1.,1.,1.,3.,2.) == Approx(-0.6454972243679027).epsilon(1e-10));
  CHECK(TACGCoe::CG(1.,1.,3.,1.,3.,2.) == Approx( 0.6454972243679027).epsilon(1e-10));
} // end TEST_CASE
#endif
